Matrix Product Density Operators: Simulation of finite-T and dissipative systems. 
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We show how to simulate numericaUy both the evolution of ID quantum systems under dissipation 
as well as in thermal equilibrium. The method applies to both finite and inhomogeneous systems 
and it is based on two ideas: (a) a representation for density operators which extends that of matrix 
product states to mixed states; (b) an algorithm to approximate the evolution (in real or imaginary 
time) of such states which is variational (and thus optimal) in nature. 
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The physical understanding of quantum many-body 
systems is hindered by the fact that the number of pa- 
rameters describing the physical states grows exponen- 
tially with the number of particles. Thus, even for a 
relatively small number of particles, most of the prob- 
lems become intractable. During the last decade, how- 
ever, several numerical methods have been put forward 
which allow us to describe certain many-body systems. 
One such method is the so-called density matrix renor- 
malization group (DMRG), which is very well suited to 
describe ID systems on a lattice with short-range inter- 
actions and at zero temperature [lj,y,y|. In fact, DMRG 
has had a very strong impact in the field of condensed 
matter physics, allowing us to describe such systems with 
an unprecedent degree of accuracy and to extract their 
physical behavior. 

Several years after its discovery, DMRG was extended 
to finite temperature ID systems 4]. The method ap- 
plies to translationally invariant infinite systems, since 
in that case the evaluation of the partition function can 
be recasted in terms of finding the maximum eigenvalue 
of a finite matrix, which in turn can be found using a 
variation of the original DMRG method to classical 2- 
dimensional spin systems [5j . The main restriction of the 
method is that it does not work well for low temperatures 
and that it cannot be applied in situations in which the 
number of particles is finite and/or not homogeneous, as 
it is e.g. the case of atoms in optical lattices |^. Re- 
cently, time-dependent versions of the DMRG method 
have been put forward 0; IS H uM ■ 

The success of the DMRG method and its extensions is 
based on the fact that the many-body states that appear 
in some ID problems can be very well described in terms 
of the so-called matrix product states (MPS) 12. 13. 14J . 
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Here, the A's are matrices whose dimension is bounded 
by some fixed number D and d is the dimension of the 
Hilbert space corresponding to the physical systems. In 
fact, DMRG can be viewed, both for finite and infinite 



dimensions, as an iterative method that for a fixed D 
determines the matrices whose state h/'mps) minimizes 
the energy in a variational sense |l2l |lj| . On the other 
hand, the method introduced by G. Vidal in the context 
of simulating the dynamics of weakly ent ang led quantum 
systems [iflj and later developed in [a, Isl lll| prescribes a 
particular way to update the matrices A as a function of 
time. This method can be extended [lOj to determine the 
evolution of mixed states by considering them as vectors 
in the d^^-dimensional space of linear operators and thus 
doubling the number of matrices A. 

In this paper, we introduce the notion of matrix prod- 
uct density operators (MPDO) which extend the MPS 
from pure to mixed states. We also present the opti- 
mal way in which the time evolution of pure and mixed 
states can be approximated within these two classes 
of states (MPS and MPDO). The corresponding evolu- 
tion in imaginary time leads to a very versatile finite-T 
DMRG algorithm, not restricted to large temperatures or 
to homogeneous systems. We conclude by showing how 
dissipative systems governed by master equations can be 
efficiently simulated. 

In [lj,ll3|! a picture was introduced to analyze MPS 
and DMRG from a quantum information perspective (see 
Fig. 2 in ^). The MPS O can be depicted as be- 
ing built up by a collection of virtual D-level systems 
paired in maximally entangled states. The MPS is ob- 
tained by applying linear maps which transform the D^- 
dimensional Hilbert space of pairs of virtual systems into 
the local Hilbert spaces associated to the physical sys- 
tems. Thus, the MPS is completely specified by these 
maps which in turn can be reexpressed in terms of the 
matrices A of Eq. (^. A central result of this paper is 
that this picture can be extended to describe mixed quan- 
tum states. The idea is as follows: instead of applying 
a linear map to the Hilbert space associated to pairs of 
virtual systems, we apply a general completely positive 
map to the corresponding operator space, which is the 
most general map allowed by quantum mechanics [l6j. 
We hence define the class of MPDO as the states which 
can be obtained by this procedure. More specifically, 
a MPDO p of iV d-level particles with Di, D2, • ' ' Dn- 



dimensional bonds is defined as 
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where Af^'"^'' are Df x Dl_^^ matrices tliat can be de- 
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Here rfj, is at most dOkD^j^i, and the matrices A^' are 
of size Dk X Dk+i- Condition © is a semidefinite con- 
straint and is sufficient to ensure that the associated map 
is completely positive. Relaxing this condition gives rise 
to operators which are not necessarily positive semidef- 
inite (i.e., the picture introduced above also allows one 
to represent general operators). Note also that for suf- 
ficiently large D, Eq. |2Jl includes any density operator 
acting on the Hilbert space. For example, the maximally 
mixed state corresponds to the choice M^'* = (Sg^^'l, and 

a pure MPS is recovered when M^'^' ^ A%® {A()* . 

It is possible to express any MPDO in terms of a (pure) 
MPS by defining the latter over a larger Hilbert space, i.e. 
by using the concept of purification [ig . To each physical 
system we associate an auxiliary system (or ancilla) with 
a Hilbert space of dimension dk^ and after choosing an 
orthonormal basis |sfe, Cfe) for these particle-ancilla pairs, 
we write the corresponding MPS state as 
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The MPDO p is obtained after tracing over the ancillas, 
i.e. p = Tra(|^)($|). It is remarkable that for many 
interesting states the ancilla's dimensions can be chosen 
to be small {dk — d) and thus the purification Q) yields 
a very efficient representation of the MPDO. Note also 
that the ancilla matrices Ak in Q) can be easily recovered 
from the matrices Mk by means of an eigenvalue decom- 
position. On the other hand, given the matrices Mk, one 
can efficiently determine expectation values as follows: 



(Oi...OAr)^ = Tr(^i,c 
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where Ek,o =Y.s 

In order to measure the error that we make when we 
approximate a density operator po by some other one 
p, one can use different quantities. For instance [l6j, 
Ep{p, po) — 1—F{p, po), where F extends the notion of fi- 
delity from pure to mixed states, it is defined as the max- 
imal overlap between all possible purifications of p and 
Po, and it is given by F{p,po) = Tr{y/y/ppoy/p). Aher- 
natively, we may simply take E^sip, Po) = Tr[(p — Po)^] 




FIG. 1: Evolution of tlie ground (pure) state of an Ising 
model after switching on a small transverse magnetic field. 
We plot the infidelity, E — 1 — \{il)(t)\'4^mps{t))\, between the 
exact solution, ip{t), and the MPS computed using our itera- 
tive method (solid) and that of Ref. [g| (dashed) , using always 
matrix product states of dimension Dk ~ 4. 



which is related to the Hilbert Schmidt scalar product in 
the space of operators. 

If we want to determine a Hamiltonian time evolution 
of a mixed state in real or imaginary time, we may in- 
stead simulate the evolution of the purification. Starting 
out from the MPS purification |^(0)) Q, we take a small 
time step and compute the next purification |'I'(Ai)) ex- 
actly. Unless the evolution is local, the dimension of the 
matrices A will increase, and we will have to find a MPS 
that approximates the exact one, as in the method in- 
troduced in 8]. However, we do this here in an optimal 
way, i.e. we compute the MPS which, for given dimen- 
sions {Dk}, maximizes the overlap with the exact state. 
This "truncated" purification is taken as our next initial 
data, and we repeat the process. At any time, the pu- 
rification can be used to reconstruct the density operator 
and to compute expectation values of observables. 

Our method is based on an iteration which resembles 
the sweeps used in standard DMRG. Let us assume that 
we have \(f>), a MPS ^ and we want to find the closest 
one, |(/)') in which the ^'s are replaced by ^'s and the 
corresponding dimensions are Dk < Dk- Starting with 
a guess of the A's (for example, the one of the previous 
step), on each step of the iteration we choose a site, k, 
and find the optimal Ak such that the overlap between 
\(f>) and the MPS is maximal. The optimal choice is given 
by the solution of the system of equations 



^a,i3 {^k)a',l3' — H^p. 

The tensors C and D are defined as 
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= Tr(A„,,^, (g, Aa,,pEk+i . . . EnEi . . . Ek-i), 



Hip = Tr(^^®A,-,0-Gfe+i...GArGi...Gfc-i), 

where Aa,p is a matrix with a single nonzero element at 
row a and column (i, and Ek = X^s ^k®i-^k)* ^ ^^"^ ^fc = 




FIG. 2: (a) Magnetization vs. transverse magnetic field of the 
thermal state of a XY model with 60 spins, for temperatures 
T = 0.05,0.5,5, and 50 (bottom to top), (b) Error in the 
density matrix of the thermal state of Eq. JUJ computed using 
MPDO, vs. temperature, for a chain with N = 8 spins, D = 8 
(solid) and D = 14, 20, 24 (dashed). 



Z^s^l ® i^k)*- Sweeping from fc = 1 to A^ and back 
several times we reach always a fixed point with maximal 
overlap |(0|0)|. The method becomes most efficient for 
open boundary conditions (Z?i = -Dat+i = 1), because 
we can then always impose that on the optimized site k 

^(ij)(ij)t = 1, Vj>fc, and (7a) 

S 

^(ij)t(i^^) = 1, Vj<fc. (7b) 

Then C^ g = 5a,a'5p^p' and the new matrix A is just 
equal to H. In order to ensure conditions (|7a|) and (|7bp 
in the following step, one can do use a singular value 
decomposition of A^. in exactly the same way as it is 
done in standard DMRG. 

We have illustrated the performance of our method for 
a simple case in Fig. ^ We chose an Ising Hamilto- 
nian, H^ =Yl,i ~ ^i^t+i +Si ^^^ ■ We begin with the 
ground state of /i = 0, and switch to ft, = 0.2 abruptly. 
We have made simulations of 6 spins with MPS of dimen- 
sion D = 4. In Fig. n we plot the error of our method 
together with that of the method introduced in |3| ■ Since 
the variational procedure is optimal at each time step, it 
achieves a better result. For the cases studied below the 
use of our method was crucial to obtain very good accu- 
racy. 

Now we apply the previous ideas to approximate the 
density operator for a system at thermal equilibrium, i.e. 
p ex e~^^ with (3 = 1/T. For /3 — (infinite tempera- 
ture), the state p is a maximally mixed state 1. For any 
other temperature we use that ^17l] 
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This means that we can obtain p by evolving according 
to the Hamiltonian in imaginary time, with time step 
At. This evolution can be performed with the optimal 
method presented above, and in this case it is sufficient 
to fix the dimension of the ancilla to d. 



Using these techniques, we have computed the equilib- 
rium states of a spin s = 1/2, XY model with transverse 
magnetic field 
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(9) 



In Fig. [2Ia) we plot the magnetization versus the mag- 
netic field for four different temperatures, A'^ = 60 spins. 
The accuracy of the matrix product representation seems 
to increase with temperature. Using smaller lattice sizes 
{N — 8 spins) and a fixed magnetic field {h = 3.0), we 
compared the exact diagonalization and the MPDO re- 
sults for different dimensions, D. As Fig. |2Ib) shows, 
no large Z?'s are needed to describe the thermal states 
accurately. 

Let us now move on to the description of the evolution 
of spin systems in the presence of decoherence. For il- 
lustrative purpose, we consider a master equation of the 
form 

|p = ^ [H„p] +lY.^fpaf ~ p) =: Lp. (10) 

i 

where H^ is the Ising Hamiltonian mentioned before and 
7 characterizes the interaction with the environment. 

Given an initial condition p(0) = po written in MPDO 
form (|2Jl, we want to find the MPO which best approx- 
imates p{t). The procedure that we propose resembles 
the one for Schrodinger equations, in that we integrate al- 
most exactly the master equation for a short time At, and 
then find the optimal MPDO of smaller dimensions which 
is closest to p{t + At). Let us denote by £(£, t) — exp(i£) 
the CPM so that p{t) = £{C,t)p{0). We wiU decompose 
this operator as follows 

£{C, At) ~ [£iCe,At/2)£{Co, At)£{Ce,At/2)f , (11) 

where £ = £o -I- £e is a splitting of the Liouvillian into 
commuting terms which act on odd, (2fc -I- l,2fc), and 
even neighbors, (2fc, 2fc + 1), respectively. 

The action of the operators £{£e,At/2) or £{Co,At) 
on a state of the form |(2Jl can be computed exactly, and 
one obtains an operator po — £{Co, At)p of the same 

form but with the substitution M^l"'"'^" M2l''^^''''^''+' ^ 

^2k''2k+i^' ^'" ^''^^ ■ Now we find the optimal operator 
PL, which has the form of Eq. (O, and is built from ma- 
trices of fixed dimension. As before, by optimal we mean 
that the distance to the original one is minimized. This 
time we will use as a measure of the distance E^is {Po, Pl) 
and thus we will not impose the condition Q (the use of 
Ep in this context is more complicated and will be ex- 
plored in a future work). This leads to an optimization 
algorithm which is very much like the one introduced be- 
fore at the level of purification, but in which now po and 
PL are treated as vectors on a d^ dimensional space. 




FIG. 3: (a) Errors in simulations of a cluster state with A^ = 
6,7 and 8 spins (solid, dash and dash-dot), under Eq. IIOL 
for 7 — 0.1, h = 0, matrices of size L — 14 and timestep At — 
0.01. (b) Decay of correlations, Czz = (crfcrf_,_i) — (erf ){crf+i), 
for the site i = N/2, on GHZ states with iV = 10, 20 and 30 
spins (dots, crosses and solid line) for L = 10. 



where d ~ d x d. This way, it is possible to use the 
same code for studying the evolution of both pure and 
mixed states. We remark that, while this method cannot 
guarantee the positivity of the truncated operator pL {t) 
obtained at the end of the algorithm, if the truncation er- 
ror E^isiPo, Pl) is kept small, it will be possible to bound 
the errors made when using pL(t) to compute expecta- 
tion values of any physical quantity. Note also that the 



matrices M, 



G C 



LkXLu 



do not need to have a size 



Lfc = D\ (which is the square of a number). 

As an illustration, we have simulated the decay of cor- 



relations of the GHZ state, \ip{Q)) 



1 

V2^ 



(|0)®^ + |1)®^) 



and of a cluster state [Ig, under the evolution dictated 
by Eq. I|1U|) . In Fig. Ol^a) we show the errors for sim- 
ulations with up to A^ = 8 spins, a size for which we 
can compute the solution numerically with direct meth- 
ods (a Crank-Nicholson which already has an error of 
the order of 10^^*^). There are two sources of errors in 
our method. One is the Trotter expansion (|ll|l . but as 
it is of order ©(At^), it can be decreased by making 
shorter time steps. The second source of error is the 
truncation of the MPDO to a lower dimension. In the 
problems that we have simulated, truncation error only 
affected the cluster state and did not grow much with 
the size of the system. In Fig. O (b) we plot some sim- 
ulations made with GHZ states of up to 50 spins, from 
which an exponential decay law for the correlation func- 
tions, c^3(i) = («+^) - (CT,f)(crf_^^) ~ e-"''*, can be 
extracted. The method is rather efficient and can easily 
handle this and bigger problems in ordinary computers. 
In summary, we have introduced the concept of ma- 
trix product density operators (MPDO) and their pu- 
rification in terms of matrix product states (MPS) in 
order to describe mixed states in quantum many-body 
systems. We have also developed a method to determine 
the MPS (MPDO) which optimally describes another one 
which is composed of matrices of higher dimensions. We 
have used those ideas to introduce an algorithm which 
allows us to study ID systems in thermal equilibrium. 



In contrast to previous methods, it applies to finite sys- 
tems, with general short range interactions (not neces- 
sarily translationally invariant) and works over the whole 
range of temperatures T g [0, oo). In fact, the algorithm 
is a very natural extension of DMRG to finite temper- 
atures because for T ^ one recovers the same MPS 
that one would obtain with standard DMRG. We have 
also shown how to simulate the evolution of a system 
in the presence of dissipation, i.e. when it evolves ac- 
cording to a master equation which involves short range 
interactions. We remark that in this last case an alter- 
native possibility is to combine the concepts introduced 
here with the "quantum jump" approach [l9j. The idea 
is to perform several realizations of the evolution start- 
ing out with randomly chosen pure states in such a way 
that their mixture reproduces the initial MPDO 0. In 
each realization the pure state is evolved according to a 
stochastic Schrodinger equation consisting of a time evo- 
lution with an effective Hamiltonian which is interrupted 
by quantum jumps jl2|. After each short time step (or 
quantum jump) the new pure state is determined using 
the method introduce here and which optimally approxi- 
mates the new state with a MPS. At the end, expectation 
values of physical observables are determined by averag- 
ing with respect to the different realizations. Which of 
these two methods is more efficient will depend on the 
number of realizations one has to perform in practice in 
the latter in order to achieve a prescribed accuracy. 
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